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SCALE-BASED IMAGE FILTERING OF MAGNETIC RESONANCE DATA 
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FIELD OF THE INVENTION 

This invention relates generally to the field of magnetic resonance imaging (MRI), 
image display, image processing and image filtering to enhance structure and reduce noise. 

BACKGROUND OF THE INVENTION 

Magnetic Resonance Imaging (MRI) has revolutionized radiological imaging of the 
internal structures of the human body. It has the advantage of being noninvasive and offers no 
known health hazards. A variety of MRI protocols are currently available, however image 
20 acquisition techniques often suffer from low signal-to-noise ratio (SNR) and/or contrast-to- 
noise ratio (CNR). Although many acquisition techniques are available to minimize these, post 
acquisition filtering is a major off-line image processing technique commonly used to improve 
theSNRandCNR. A major drawback of filtering is that it often diffiises/blurs important 

Structures along with noise. 

Two- and higher-dimensional images are currently available through sensing devices 
that operate on a wide range of frequency in the electromagnetic spectrum - from ultrasound to 
visible hght to X- and y-rays (Chu et al, Foundation of Medical hnaging, John Wiley Sons, 
Inc New York, NY, 1993). The usefulness of any of these modalities depends on the 
perceptability of certain features in the created images. Acquired images are often degraded by 
30 various types of artifacts resulting in a lowering of signal-to-noise ratio (SNR) and contrast-to- 
noise ratio (CNR). Despite high adaptivity and robustness of the human visual perceptual 
system, low SNR and CNR often degrade the information contained in the image and its utility. 
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Further, image artifact affect many image processing tasks such as segmentation (Sahoo et al. 
Computer Vision Graphics and Image Processing 41:233.260 (mS); Udupae^ a/., Graphical 
Models and Image Processing 58(3):246-261 (1996)). registration (Wells et al.. Medical Image 
Analysis 1:35-51 (1996); hMrAzet al. Medical Image Analysis 1:151-161 (1996)). and 
5 visual rendition (Udupa, Radiographics 19:783-806 (1999); Hohne al, (eds.), 3D Imaging in 
Medicine: Algorithm^- Systems. AppHcations . Springer-Verlag, Berlin, 1990), which are 
crucial in many applications. Therefore, noise reduction, that is, the improvement of SNR and 
CNR, is of considerable interest in many imaging applications. 

' There are basically two approaches for improving SNR and CNR: by changing the 
10 method of image acquisition and by utihzing post-acquisition image processing. Improving the 
method of acquisition usually entails lengthening the acquisition time, or sacrificing spatial 
resolution, or upgrading of the imaging device hardware. Filtering is an off-line image 
processing approach, which, if properly designed, may be less expensive than, and as effective 
as, acquisition improvement approaches without having to sacrifice spatial resolution. 
15 ' The aim of filtering is to enhance wanted (structure) information and to suppress 
unwanted (noise) information. Filtering operations may be classified into two categories - 
enhancing (also know as high pass filtering), wherein wanted information is enhanced 
hopefiiUy without affecting unwanted information, and suppressing (also known as low-pass 
filtering), wherein unwanted information is suppressed hopefiiUy without affecting wanted 
20 information. 

Noise reducing filtering operations may be classified into two major categories. The 
first is space-invariant filtering (Rosenfeld and Kak. Digital Picture Processing. Academic 
Press. Inc.. Oriando. FL, 1982). wherein a spatially independent fixed smoothing operation is 
perfomied over the entire image. The other is space-variant filtering (Lee, IEEE Transactions 
25 on Pattern Analysis and Machine Analysis 2:165-168 (1980)), wherein the smoothing operation 

is modified by local image features. 

A major drawback of space-invariant filtering techniques is that, along with noise, they 
often blur important structures. To overcome this problem, a variety of local image feature- 
dependent adaptive filtering strategies have been developed, including local image-statistics- 
30 driven filtering (Lee. 1980), least-squares error filtering (Chan and Lim, IEEE Transactions on 
Acoustic, Speech, and Signal Processing 33:1 17-126 (1985)), gradient inverse-weighted 
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filtering (Wang et al.. Computer Graphics and Image Processing 15:167-181 (1981); Wang, 
IEEE Transactions on Signal Processing 40:482-484 (1992)). multiple model Kalman filtering 
(Woods et al., IEEE Transactions on Pattern Analysis and Machine Analysis 9:245-253 
(1987)), adaptive recursive filtering using local image features (Rank et al., IEEE Transactions 
5 on Image Processing 1:431-436 (1992)), local shape-based template-matched adaptive filtenng 
(Ahn al., IEEE Transactions on Medical Imaging 18:549-556 (1999)), and gradient- 
controlled anisotropic diffusive filtering (Perona et al., IEEE Trans-actions on Pattern Analysis 
and Machine Analysis, 12:629-639 (1990)). This later method has become quite popular, and 
has been fiirther extended and utilized (Gerig et al. , IEEE Transactions on Medical Imaging 
10 1 1 :221-232 (1992); Jackson et al, J. Computer Assisted Tomography 17:200-205 (1993); Bajla 
et al, Proc. MEDINFO, Part 1:683-686 (1995); Parker et al, Proc. International Society of 
Magnetic Resonance in Medicine 7:175(1999)) in a variety of applications. 

Anisotropic diffixsive filtering (Perona et al, 1990) provides a method, in which m each 
iteration, intensity diffiision takes place between every two spatially close image spatial 
15 elements ('spels'). THe magnitude of the diffusion between the two spels is determined by the 
magnitude of the gradient of image intensity between them. It varies in a monotonically, non- 
decreasing manner up to a certain gradient magnitude value, and subsequently changes m a 
monotonically, non-increasing fashion. The flow across boundaries is thus significantly 
reduced by utilizing information about the presence of boundaries provided by intensity 
20 gradients. 

A significant drawback of this very usefiil approach, however, is that it does not offer 
any image-dependent guidance for selecting the optimum gradient magnitude. More 
importantly, since it does not use any morphological or structural information to control the 
extent of diffixsion in different regions, fine structures often disappear and fixzzy boundanes are 

2 5 further blurred upon filtering. 

Thus, until the present invention, there has remained a need in the art for a reUable 
method for filtering SNR and CNR in an MRI acquired image in which intensity gradients are 
accurately determined, but in which selection of the optimum gradient magnitude is permitted 
and fine structures are detected. To overcome the problems of diffusive filtering in the pnor 

30 art, the inventors have explicitly utilized 'object size' information - or 'scale' - to control the 
degree of smoothing that is done in different regions of the image. The notion of scale, as 
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defined herein, is different from the term 'scale' used in scale-space computer vision literature 
(Lindeberg, Scale-Space Theory in Computer Vision . Boston, MA: Kluwer. 1994; Lifshitz et 
al, IEEE Transactions on Pattern Analysis and Machine Intelligence 12(6):529-5407 (1990); 
McAuliffe et al., in Proc. Visualization in Biomedical Computing, Lecture Notes in Computer 
Science, 1 131, Springer Verlag, Berlin, 1996, pp. 173-182; Pizer et al.. Computer Vision and 
Image Understanding, 69(1):55-71 (1998)), although the intent is the same, i.e., to tailor the 
processing to local object size. 



10 



15 



20 



25 



30 



SUMMARY OF THE INVENTION 

The invention provides novel space-variant anisotropic (scale-based) methods for 
filtering SNR and CNR in an MRI acquired image - a spatial-resolution adaptive scale- 
computation method, and a scale-based neighborhood averaging method, and a scale-based 
di ffixsive^filtj ring method -both of the scale-based methods outperform the anisotropic 
diffusive approach in the prior art. The present scale-based filtering methods use local 
structure size or "object scale" information to arrest smoothing around fine structures and 
across even low-gradient boundaries. One method uses a weighted average over a scale- 
dependent neighborhood, while another employs scale-dependent diffiision conductance to 
perform filtering. 

Both methods adaptively modify the degree of filtering at any image location, 
depending on local object scale. This permits the accurate use of a restricted homogeneity 
parameter for filtering in regions with fine details and in the vicinity of boundaries, while at the 
same time, a generous filtering parameter is used in the interiors of homogeneous regions. 

Scale-based filtering basically uses the local scale information at every spel to control 
the extent of filtering. By definition, the scale in regions with fine details or in the vicinity of 
boundaries is small. Thus, a restricted homogeneity parameter is used for filtering in small- 
scale regions (corresponding to fine structures and vicinity of boundaries), and to use a 
generous filtering parameter at large-scale regions (corresponding to interiors of large 

homogeneous regions). 

Consequently, the local scale information indicates the radius of the largest hyperball of 
a "homogeneous" region centered at the spel. This concept was originally introduced by the 
inventors at the SPIE Conference Medical Imaging 2000, and fiirther described in Computer 
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Vision and Image Understanding, 77:145-174 (2000), where it was objectively demonstrated 
that the method significantly improve image segmentation based on fuzzy connectedness. 
Thus, the present invention provides a considerable body of evidence to show the effectiveness 
of the scale information in improving diffusive filtering, as well as for devising simple new 
5 strategies based on weighted averaging of intensities. 

Qualitative comparison experiments based on both mathematical phantoms and medical 
(patient) MR images utilizing the anisotropic diffusion method and the two scale-based 
methods, showed significant improvements using the scale-based methods over the extant 
anisotropic diffusive filtering method in preserving fine details and sharpness of object 
10 boundaries. Quantitative analysis on geometric phantoms generated under a range of 

conditions of blurring, noise, and background variation confirmed the superiority of scale- 
based approaches embodied in the invention. 

Thus, a scale-based filtering method of the invention is provided, wherein an enhanced 
MRI-acquired image is achieved for any selected MRI protocol, and for any selected region of 
a patient's body, by the steps comprising: imaging the region of the patient's body by the 
selected MR protocol to form an acquired image; and filtering the acquired image by a scale- 
based spatial resolution adaptive method using region-constancy based on local homogeneity to 
produce an enhanced image. The resulting image is enhanced in a manner independent of 
variations within or between patients, within or between tissues being imaged, or within or 
between MR devices used to acquire the image. 

Also provided are computerized systems and devices for implementing the foregoing 
methods. 

In addition, as provided in at least one embodiment, scale-based filtering methods of the 
invention are automated. Moreover, an embodiment of the invention is provided in which 
scale-based filtering methods are built into an MR scanner, permitting production of enhanced 
real time images. 

Also provided is an enhanced MR image for imaging a region of the body of a patient 
according to one or more of the embodied scale-based filtering methods of the invention, 
including a scale-based neighborhood averaging method, or a scale-based diffusive filtering 
method. 
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under a variety of conditions. invention will be set forth in 

invention. 
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flow magnitude function that complies with the principle of nonlinear diffusion. FIG. 3(c) 
graphically depicts the shape of different conductance functions at different scales. The right- 
most curve of FIG. 3(c) represents the conductance function at maximum scale tmax- FIG. 3(d) 
graphically depicts the shape of different diffusion flow magnitude functions at different scales. 
5 The right-most curve of FIG. 3(d) represents the diffusion flow magnitude at maximum scale 

^MAX. 

FIGs. 4(a)-4(h) depict comparisons among the three filtering methods using the 2D slice 
of a MR 3D scene of the head of a multiple sclerosis patient originally shown in FIG. 1(a). 
FIGs. 4(a) and 4(b) depict 2D scenes within two ROIs, one fi-om the top left comer and the 
10 other fi-om the bottom right comer of the original scene in FIG. 1(a). FIGs. 4(c)-4(h) depict 
filtered scenes (in 2D) corresponding to the scenes in FIGs. 4(a) and 4(b) that resulted from 
anisotropic diffusion (shown in FIGs. 4(c) and 4(d)), from scale-based neighborhood averaging 
(shovm in FIGs. 4(e) and 4(f)), and fi-om scale-based diffusion (shovm in FIGs. 4(g) and 4(h)). 



1 5 scene of the head of a multiple sclerosis patient. FIG. 5(a) depicts a 2D slice of the original 3D 
MR scene and the selected ROI. FIG. 5(b) depicts the scene within the ROI, magnified. FIGs. 
5(c)-5(e) depict filtered scenes corresponding to the scene in FIG. 5(b) that resulted fi-om 3D 
filtering using anisotropic diffusion (shown in FIG. 5(c)), using scale-based neighborhood 
averaging (shown in FIG. 5(d)), and using scale-based diffusion (shown in FIG. 5(e)). 

20 FIGs. 6(a)-6(f) depict comparisons among the three filtering methods using a 2D 

phantom scene containing different geometric shapes of different size and orientation. FIG. 
6(a) depicts an original phantom scene (one among 25 such scenes) at a low level of blurring 
and a medium level of noise, together with an intensity profile along the depicted bright 
horizontal line. Six ROIs are also shown within FIG. 6(a), which are used for FIG. 7. FIG. 

25 6(b) depicts a magnified version of a section of the profile in FIG. 6(a), which section is shown 
at the bottom of FIG. 6(a) as a rectangle outlined in black. FIG. 6(c) depicts the profile 
corresponding to that in FIG. 6(b), but which were obtained for the original binary phantom 
scene after adding only the background ramp component (but not adding blurring and noise). 
FIG. 6(d)-6(f) depict profiles corresponding to that in FIG. 6(b), but which were obtained from 

30 the filtered scenes: by anisotropic diffusion (shown in FIG. 6(d)), by scale-based neighborhood 
averaging (shovm in FIG. 6(e)), and by scale-based diffusion (shown in FIG. 6(f)). 



FIGs. 5(a)-5(e) depict comparisons among the three filtering methods using a 3D MR 
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FIGs. 7(a)-7(d) depict comparisons among the three filtering methods using the 2D 
phantom scene of FIG. 6(a) for the six ROIs shown therein. FIG. 7(a) depicts scenes within the 
six ROIs taken from the unfiltered phantom scene. FIGs. 7(b)-7(d) depict the corresponding 
fihered scenes that resulted from anisotropic diffusion ((shown in FIG. 7(b)), from scale-based 
5 neighborhood averaging (shown in FIG. 7(c)), and from scale-based diffusion (shown in FIG. 
7(d)). 



DESCRIPTION OF PREFERRED EMBODIMENTS OF THE INVENTION 

The invention provides novel scale-based MR image filtering methods, including scale- 
10 based neighborhood averaging and scale-based diffusion, which avoid or prevent blurring 
around fine structures and across boundaries using local structure size or "scale" information. 
Fundamentally, local structure size, or "scale", as defined by region homogeneity, improves the 
filtering operation by its utility in controlling the operation. The scale concept allows a 
jl restricted smoothing selectively to be achieved in the vicinity of boundaries and in fine 

(, ; 

Ci 1 5 structured regions, while at the same time allowing more generous filtering in the interiors. 

Thu s, ^^cale^b^ed fi ltering method of the invention is provided, wherein anenhanced ^^^1-^ 
hj acquired image is achieved for any selected MRI protocol, and for any selected region of a 

- patient's body, by the steps comprising: imaging the region of the patient's body by the 

C3 selected MR protocol to form an acquired image; and filtering the acquired image by a scale^ 

20 based spatial resolution adaptive method using region-constancy based on local homogeneity to 
2f produce an enhanced image. 

hi In the discussion which follows, both scale computation and the filtering operation are 

set up in a general n-dimensional framework, which also takes into account the anisotropy of 
spatial resolution. The framework is designed in such a way that there is only one parameter to 

25 be specified - the sole image-dependent homogeneity parameter - and methods are provided 
for selecting this parameter utilizing either user guidance or global image statistics. However, 
to better imderstand the present disclosure, certain terms and notations are defined as follows. 

Most digital imaging systems produce images on a rectangular grid where the basic 
image elements, referred to as "spels," are hypercuboids. The term "spels" is an abbreviation 

30 for spatial elements; in 2D spels are 'pixels,' while in 3D they are ^voxels.' Mathematically, 
this is explained by letting n-dimensional Euclidean space be divided into hyper cuboids by 
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n Orthogonal families, each of equally distant parallel planes. Let vbe an w-tuple whose / 
element v, represents the distance between two successive parallel planes perpendicular to the 
i^^ coordinate axis. (A 'tuple' is a recognized term, simply referring to a mathematical relation 
between its attributes.) In other words, v, is the size of each spel along the direction of the 
5 coordinate axis. 

The term v as used herein means the resolution of the grid system. Further, a coordinate 
system is chosen such that the centroid of each spel has the coordinates of the form (c\ vi, C2 vi, . 
' ' yCn Vn) where c/'s are integers. Thus, a one-to-one and onto mapping between the set of spels 
and as follows: an /i-tuple c, denoting a point in Z", is mapped to a spel, such that cj vj, for 1 
10 <j<n, gives the coordinate of the center of the spel represented by c. With this 

interpretation of spels, itself represents the set of all spels of 9?'', and concepts of spels and 
points denoting their centroids may, therefore, be used interchangeably to describe the 
invention. 



15 experiments; however, fuzzy relations for a should also be considered, especially those 
functional forms for which closely resemble the point-spread function of the imaging 
device. A fuzzy relation or on Z" is said to be a "fuzzy spel adjacency" if it is reflexive and 
symmetric (for a description of 'fuzzy relations,' see Kaufinann, Introduction to the Theory of 
Fuzzy Subsets , 1, Academic Press, New York, 1975). It is desirable that a be such that its 



20 membership value //^(c, d), for any c, J e 2", is a non-increasing function of the distance \\c-d\\ 
ri between c and d, where || • || denotes any L2 norm in 91 . 

The functional form of that is used for the image processing results presented herein, 
is as follows (however, a general functional form for jUa is considered for all formulations). For 
any spel c, /udc, c) = 1. Further, for any spels c and d, jij^c, d)^\,\ic and d differ in exactly 
25 one coordinate by 1. Otherwise, j^dc, d) = 0. 



For computational simplicity, a hard relation of adjacency was used for a in all 



In 2D: 




(Equation 1) 



30 
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Thus, it is readily seen that a is indeed a fuzzy spel adjacency. The triplet (Z", a, v) is, 
therefore, called a "fiizzy digital space," in which a is any fuzzy spel adjacency, and vis the 
resolution of the grid system. A fuzzy digital space is a concept that characterizes the 
underlying digital grid system, independently of any image-related notions. 

5 A "scene over a fuzzy digital space" (Z", a, v) is a pair C = (C,/), where C= {c | -bj 

< Cj < bj for some beZ"}, Z" is the set of w-tuples of positive integers, /is a function whose 
domain is C, called the "scene domain", and whose range is a set of integers [Z, H]. The terms 
"scene" and "image" are used interchangeably, herein, as recognized in the art. The invention 
generalizes in a straightforward manner to vector- valued scenes, z.e., when /is a vector- valued 

10 function. However, for the sake of simpUcity of description, /is treated herein as a scalar- 
valued function. A scene ^over the triple (Z", v) is referred to as a "binary scene" if the 
range of/is {0, 1}. 
1) Anisotropic Diffusive Filtering. 

A diffusive filtering process can be defined using the divergence operator 'div' on a 

1 5 vector field. A mathematical formulation of the diffusion process on a vector field V at a point 
c in coordinate-fi*ee form is given by: 

— = div V = lim f V • ds, (Equation 3) 

dt ^ 

20 wherein Ar is the volume enclosed by the surface s (surrounding c), and ds = uds, wherein u 
is the unit-vector normal and outward-directed to the 'infinitesimal' surface element ds. The 
diffusion process is controlled by the intensity flow vector field V. The functional form for 
flow is defined as 

25 V = GF, (Equation 4) 

wherein G is the diffusion-conductance function, and F is the scene intensity gradient vector 
field. 

In an isotropic diffusion, the conductance function G is constant, and it has been argued 
30 by Perona et al, 1990, that such diffusion methods blur object boundaries and structures. As a 
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result, Perona et al, proposed an anisotropic diffusion method in which G is spatially varied as 
a nonlinear function of the magnitude of the scene intensity gradient, such that the conductance 
function encourages smoothing within a region (characterized by low intensity gradients) in 
preference to smoothing across boundaries (characterized by high intensity gradients). 
5 2) Scale-Based Filtering. 

Although, the main motivation for anisotropic diffusive filtering was to reduce noise 
while minimizing blurring across boundaries, since it does not account for local structure size, 
it often loses/blurs fine structures and diffuses across boundaries. This phenomenon is 
illustrated in FIGs. 2, 4-7, which will be discussed in greater detail below. 

10 In this section, however, two scale-based filtering methods are presented that use the 

structure size information to accurately arrest diffusion around fine structures, and across even 
low gradient boundaries. The size of local structures is determined at every spel under a pre- 
specified scene-dependent, region-homogeneity criterion, and then the generosity of filtering is 
controlled using this parameter. For example, in FIG. 1(a), which is a 2D slice from a MR 3D 

15 scene of the head of a multiple sclerosis patient, using a suitable region-ho|mogeneity 

parameter, a local determination is done to determine the largest disc that can be centered at 
any point, such as /?, within which the intensities are homogeneous. In FIG. 1(a), the size of the 
local structure to which pixel p belongs is bigger than that to which q or o belong. The scale is 
usually small in the vicinity of boundaries, such as at o. It is, therefore, reasonable to allow a 

20 more generous filtering at p (/.e., an averaging over a larger neighborhood, or a diffusion using 
a higher conductance), than that at q or at o. 

The theory of defining scale as the size of the local structure and utilizing scale in 
carrying out local operations was originally proposed by the inventors (Saha et al, 2000, herein 
incorporated by reference). Following the same line of reasoning, "scale" in a scene Cat any 

25 spel, cgC, is defined as the radius, r(c), of the largest hyperball centered at c, which lies 

entirely in the same object region (defined under a pre-specified region-homogeneity criterion). 
Ironically, it appears that object definition is first needed to determine scale. Thus, a simple 
and effective algorithm (presented in the following section, as Algorithm OSE) has been 
defined by the inventors (see also Saha et al, 2000), which estimates r{c) at every spel, ceC, in 

30 any scene C without explicit object segmentation, but based on continuity of region 

homogeneity. This 'rough' object size information is utilized to control the diffusion process 
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SO as to improve diffusive filtering and to devise a simple neighborhood averaging strategy that 
can also out-perform diffusive filtering. 
2.1) Computation of Scale. 

A "hyper ball Bk,iic) of radius k>0 and v^ith center at c e C in a scene C= (C, J) over 
(Z^ a, v) is defined by: 



(Equation 5) 



wherein the term e is defined within the equation as a spel belonging to C 

Note that 5k,v/(c) forms a digital hyper ball when the spatial resolution is isotropic, and a 
10 digital hyper ellipsoid when it is not. The equation is formulated so that, in both cases, 5k,v/(c) 
represents a (quantized) hyper ball in the scene domain. For the purposes of the invention, 
Bk,Xc) is referred to as a hyper ball, irrespective of whether or not vis isotropic. 

For a hyper ball Bk,xic) of any radius A>0, and centered at c, a fi-action, FOk,yic) 
("fi-action of object"), is defined to indicate the fi*action of the ball boundary occupied by a 
1 5 region which is sufficiently homogeneous with c, by: 



|s.,vW-fi.-,..W| 



wherein 1 Bk,yic) - BkA.xip) \ is the number of spels in {i.e., the volume of) Bk,Jic) - BkA.J^c), 
20 and Wy^ is a homogeneity fimction. A detailed description of the characteristics of is foimd 
in Saha et al, 2000, herein incorporated by reference. The membership fimctions utilized in 
fiizzy subset theory corresponding to the fiizzy proposition "x is small" can be used for Wy;. A 
zero-mean, unnormalized, Gaussian fiinction is a preferred choice for Wy^ 

The algorithm for object scale estimation is presented below as Algorithm OSE, 
25 wherein "OSE" refers to "Object Scale Estimation." It has been modified firom that which was 
originally presented by Saha et al, 2000, to account for resolution anisotropy of acquired 
scenes. 

Algorithm OSE 
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Input: C,c e Cy a fixed threshold ts. 

Output: r(c). 

begin 

SQtk= 1; 

5 while FOk^Jic)>ts do 

set A: to A: + 1 ; 
endwhile\ 
set r{c) to A: - 1; 
output r(c); 
10 end. 

The algorithm iteratively increases the ball radius A: by 1, starting fi:om 1, and checks for 
FOk,iic), the fi"action of the object containing c that is contained in the ball. The first time that 
this fi"action falls below tsy the ball then contains an object region different fi"om that to which c 
belongs. 

15 In accordance with Saha et al, 2000, ts = 0.85 is used. The rationale for this choice is 

as follows: In a 3x3 neighborhood of a pixel c, in a 2D scene, Q one out of the eight 
neighboring pixels may belong to a different object (to account for noise), but the 
neighborhood is still considered to be entirely within the same object. This leads to ts = %. No 
significant change in the results was observed by varying ts in the range 0.8 < < 0.9. 

20 To complete the description of scale computation, the following is a description of two 

methods by which a proper value is selected for region-homogeneity parameter cr^^r. In the first, 
which is an operator interactive training-based method, an operator using a mouse controlled 
brush, paints on a display of a slice of an image to specify objects of interest, excluding inter- 
object boundaries. The mean, Af/,, and the standard deviation, an, of local intensity gradient 

25 magnitudes, are computed as the mean and the standard deviation of the intensity differences, 
\f{c) -fid) I , for all possible pairs (c, d) of spels in the painted regions, such that c^d and 
fiaicyd) > 0. These estimates are then used to set up the value of the homogeneity parameter as 
follows: 

30 CT^ = Mh-^thah, (Equation 7) 
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wherein th is a constant, and its value is set to 3. This is based on the rationale that in a normal 
distribution, 3 standard deviations on the two sides of the mean cover 99.7% of the population. 

The second method of selecting is similar to the preceding method, except that it 
requires no operator intervention. Instead of computing the mean and the standard deviation of 
5 local intensity gradients over the painted regions, the local, in-plane intensity gradients (/.e., c 
and d are in the same plane) are first determined over the entire scene domain. Local, out-of- 
plane intensity gradients are excluded because, in most medical imaging applications, the out- 
of-plane resolution is relatively coarse. However, local out-of-plane intensity gradients may be 
used when the spatial resolution is isotropic. The upper 10th percentile of the gradients are 

10 then discarded in order to accovmt for inter-object.boundaries. The mean and the standard 

deviation are computed over the lower 90 percentile of the gradients. Finally, the expression in 
Equation 7 is used for determining <jy/. In this fashion, the only parameter required in the 
method is determined automatically. 

In the examples that follow, the second method was used for computing a^^. For 

15 example, the scale image shown in FIG. 1(b) corresponds to the MRI slice of the head of the 
patient shown in FIG. 1(a). In the scene in FIG. 1(b), intensity at a location is proportional to 
the scale value at that location. As can be seen, the interior regions have higher scale (brighter 
intensities), while in contrast, the detailed structural regions, as well as the regions in the 
vicinity of boundaries, have lower scale (darker intensities). 

20 2.2) Scale-Based Neighborhood Averaging. 

A fixed (Gaussian) kernel-based filtering method essentially convolves the original 
scene with a fixed kernel. A major drawback of such a method is that, when the kernel is 
defined over a small neighborhood, the method fails to significantly increase the SNR upon 
filtering. On the other hand, when the kernel is defined over a large neighborhood, smaller 

25 shapes are lost and substantial amount of blurring takes place across boundaries. Finding an 
optimum neighborhood size (kemel) had been a fiindamental problem in this type of approach 
in the prior art. By bringing object scale into the fi-amework, not only is the size of the 
neighborhood specified, but also it is allowed to vary over the scene domain in a meaningful 
way. This allows large neighborhoods to be used in the interiors of large structures (e.g., at p in 

30 FIG. 1(a)), while small neighborhoods are used in fine, detailed regions and in the vicinity of 
boundaries (eg., at q or at 6) in the same figure. 
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A scale-based neighborhood averaging method can be formally defined as creating an 
output filtered scene C^, = (C,/sa.) for an input scene C= (C,/) over (Zny a, v), where, for any 



wherein coc is a distance-dependent weighting function. The true distance between c and e is 
used for the weighting function cocy with the purpose being to develop a scale-dependent 
filtering kernel using Br{c)A^)* The filtering kernel cDc should satisfy the following properties: 

1 . The domain of cd^ is the set of non-negative numbers. 

2. The range of o)cy is [0,1] and o)c{0) = 1. 

3. There is one control parameter a^^ of cDc, such that for a given distance between c and 
e, cOc, is monotonically non decreasing with <j^^ . 

4. For any given control parameter , cOc is monotonically non-increasing with the 
distance between c and 

In the fuzzy subset theory, the membership functions corresponding to the fuzzy 
proposition "x is small" satisfy these properties. Three such popular functions are (1) 
unnormalized Gaussian, (2) ramp, and (3) box window. In the examples that follow, an 
unnormalized Gaussian function for fit was used. Intuitively, it follows that cr^^ should be a 

function of r(c), and accordingly = r{c) is used. 
2.3) Scale-Based Diffusive Filtering. 

In the original anisotropic diffusion method, flow across objects is reduced due to high 
intensity gradients at boundaries. Intensity gradients at spels on blurred boundaries, where 
gradient is distributed over a thick boundary band, are often not high enough to arrest flow. As 
a result, diffusion flow takes place across blurred boundaries. The effect of this phenomenon 
becomes visually prominent when the filtering method is executed for several iterations. 
Within an 'iteration,' estimated intensity flow diffuses between every two adjacent spel in the 
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image domain. This effect is demonstrated in FIG. 2(d), wherein the diffusion process is run 
for 10 iterations on the 2D scene shown in Figure 2(b). 

Thin and fine structures in images often fade or disappear, even after fewer than 10 
iterations of anisotropic diffusion. The behavior of anisotropic diffusion around thin and fine 
5 structured regions is fiirther discussed in the Examples that follow; see also FIGs. 4, 5 and 7. 

The purpose behind scale-based diffiision is to effect a restricted diffusion around thin 
and fine structured regions and in the vicinity of boundaries, while at the same time allowing a 
generous diffusion in the interiors of homogeneous regions. This is achieved in the present 
embodiments by using a variable diffusion conductance fiinction that is controlled by the scale 

10 parameter. As demonstrated in the section entitled Scale Based Filtering and shown in FIG. 1, 
scales in the vicinity of boundaries and in thin and fine structured regions are small, while with 
the interiors they are large, thus showing the effectiveness of a scale-dependent diffiision 
process. Compare FIG. 2(c) with 2(e), which were obtained after 3 iterations, and FIG. 2(d) 
with 2(f), which were obtained after 10 iterations. See also FIGs. 4, 5, and 7. 

15 However, before fiirther describing the scale-based method, a discussion of the 

characteristics of diffiision conductance fimctions is needed. 

Perona et al, 1990, described qualitatively the shapes of the diffiision conductance 
fimction, G, and diffiision flow magnitude fimction, |V| (FIGs. 4 and 6 in Perona et al, 1990) 
that comply with the idea behind nonlinear diffiision. The required shapes of diffiision 

20 conductance and flow magnitude fimctions are shown in FIGs. 3(a) and (b), respectively. The 
necessary properties of G may be formulated as follows: 

1 . The domain of G is the set of all non-negative numbers. 

2. The range of G is [0, 1], and G = 1 when |F| = 0, while G -> 0 when |F| oo. 



monotonically increasing for |F| < cr, while it is monotonically decreasing for |F| >a. 
Intuitively, it appears that G should be monotonically decreasing with |F|. However, not 
every monotonically decreasing fimction results in a valid flow fimction. For example, if G is 



25 



3. 



There is a control parameter cr, such that for any given |F| and a\ < 02, G^^ <G^^ . 
Further, |V| (see Equation 4) has a maximum value at |F| = cr, and |V| is 



selected in the expression G = 



1 



the resulting conductance fimction is monotonically 
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decreasing, but it does not result in a valid flow function because the latter is a throughout- 
increasing function of |F|, and it possesses no finite optimum. Therefore, care must be taken in 
selecting a conductance function so that the resulting flow function possesses a unique 
maximum. Two valid conductance functions were proposed in Perona et al, 1990, and Gerig et 
ai, 1992, 

^i(F) = ^^Tm I>^>0' (Equation 9) 



1+ 



_jaL 

and G2 (F) = e^^\ (Equation 1 0) 



Both functional forms in both anisotropic and scale-based diffusion methods have been 
applied to patient MR images, and no significant visual differences were produced in the 
results. Consequently, the unnormalized Gaussian function of Equation 10 was used herein for 
all evaluations. 

1 5 As described above, cris the control parameter that determines the generosity of 

filtering. When cris large, the generosity of filtering is high, and the possibilities of blurring 
across boundaries increase. On the other hand, when cris small, the generosity of filtering is 
low, and more noise survives after filtering. Following this principle, a may also be thought of 
as a region-homogeneity parameter. 

20 As previously pointed out, the basic premise behind scale-based diffusion is to allow 

only a restricted diffusion at small-scale regions, while at the same time being generous at high- 
scale regions. This is achieved by using an adaptive homogeneity parameter cr^, in place of crin 
Equation 10, as a monotonically, non-decreasing function of local scale. Consequently, oi was 
used herein as a linear function of scale. 

25 Given this understanding of scale-based diffusion conductance, the process of scale- 

based diffusion is completely described by Equations 3, 4, and 10. Accordingly, throughout the 
remainder of the disclosure, the details of the scale-based diffusion process are presented as 
applicable to a scene over (Z„, a, v). Like anisotropic diffusion, scale-based diffusion is also an 
iterative process, and let / denote the iteration number. Let C= (QJ) be the given scene, and 
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let Q = (C,yj) denote the scene resulting from the diffusion process at the iteration. For any 
c,d e C, such that c 9^ and /j^c, cO ^ 0. let D(c, d) denote the unit vector along the direction 
from the point c toward the point d. Let F,(c, d) represent the intensity gradient vector along 
D(c, d) given by 

¥Xc,d) = -p^M^^ffi=D(c,rf) (Equation 1 1) 



min 



Note that the preceding expression takes into account resolution anisotropy. 

While determining the intensity flow from spel d to spel c, the effective scale, reff{c, d), 
10 defined as follows, was utilized. 

reffic, d) = min r{d), r{e)] (Equation 12) 

wherein e e C is the neighboring spel of c, just opposite to d, defined by e = c - (c? - c). When 
15 c is in the border of the scene domain, e was let equal to c. Using this definition of effective 
scale, the diffusion conductance function, Gt{c, d) for the flow from dtoc^t the t^^ iteration, is 
defined by: 



GXc,d) = e (Equation 13) 

wherein oi(c, d) is the scale-adaptive region-homogeneity parameter used for controlling the 
intensity flow from d to c, given by 



aXc.d) = ^I1±JL^ (Equation 14) 



In this case r^Ax is an upper limit on the scale in Q and in the following examples r^Ax ^ 5. 
(Jif, is the homogeneity parameter defined above in Section 2.1, entitled Computation of Scale . 
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Accordingly, intensity flow vector V/ (c, d) from spel c to spel d at the iteration is 
defined by: 

V,(c, d) - F,(c, ^ G/ (c, d) (Equation 1 5) 

5 Note that V, (c, d) D(Cy d) is positive when ft{c)>ft {d)\ it is negative otherwise; and it is zero 
when c = d or ft{c) = ft(d). Thus, the flow direction between two spels is always such that it 
tries to reduce the gradient between them. FIG. 3(c) shows diffusion conductance functions at 
different scales while FIG. 3(d) shows the magnitude of intensity flow vector at different 
scales. As shown in the two figures, at the maximum scale r^xx, conductance is highest and the 
10 flow magnitude is maximum. Moreover, both diffusion conductance and flow magnitude 
decrease as scale decreases. 

With these formulations, the scale-based iterative process is defined by: 



15 



20 



r/(c), for^ = 0, 

^*'> = U,-^.E..c/'.fe'')V,-,fe<*) D(c,rf). <>0. (Equation 16) 

where Ko is a constant. It can be shown in a manner analogous to the proof in Gerig et aL, 
1992, that in order to guarantee a "monotonic variation of spel intensities" (z.e., the intensity of 
any spel ceC is either non-increasing or non-decreasing with iteration t, but is not oscillatory), 
the 'upper limit' of Ko should satisfy the following inequality. For any ce Z^, 



^ ^ ^ , (Equation 17) 

Although Kd has no lower limit, more iterations are needed to achieve the same degree 
of filtering when a small value for is used for Ko^ In the following examples, the value Ko was 
25 used at its upper limit. Thus, according to the adjacency criterion (see Equations 1 and 2), Ko = 
V5 in 2D mdKo=^hin 3D. 

The embodied methods of the invention are useful when implemented by recognized 
techniques to produce an automated, scale-based filtered image, and are particularly useful 
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when installed by recognized methods in the art into an MR scanner, permitting production of 
enhanced real time images. 

Also provided in the invention is an enhanced MR image for imaging a region of the 
body of a patient according to one or more of the embodied scale-based filtering methods of the 
invention, including a scale-based neighborhood averaging method, or a scale-based difftisive 
filtering method. Further provided in the invention are computerized apparatus and devices 
incorporating one or more of the embodiments of the invention, including one or more of the 
embodied scale-based filtering methods of the invention to produce an MR image in which fine 
details and sharpness of object boundaries are preserved under a variety of conditions. 

EXAMPLES 

The present invention is fiirther described in the following examples in which both 
qualitative and quantitative experiments were conducted to assess the relative performance 
5 among the three methods - (1) anisotropic diffusion, (2) scale-based neighborhood averaging, 
and (3) scale-based diffusion. These examples are provided for purposes of illustration to those 
skilled in the art, and are not intended to be limiting unless otherwise specified. Moreover, 
these examples are not to be construed as limiting the scope of the appended claims. Thus, the 
invention should in no way be construed as being limited to the following examples, but rather, 
10 should be construed to encompass any and all variations which become evident as a result of 
the teaching provided herein. 

The principles and algorithms for both scale-based neighborhood averaging and scale- 
based diffusion were as described above. Also as described, were how to determine the sole 
image-dependent parameter, Gxf/, and how to fix other image-independent terms. For the 
15 anisotropic diffiision method, essentially the method of Gerig et aL, 1992 was followed, except 
that in computing the image gradients, resolution anisotropy was considered as in Equation 1 1 . 

All experimental results presented in this paper were obtained using the functional form 
for diffusion conductance of scale-based neighborhood averaging in Gerig et aL, 1992, which 
was the same as Equation 10, above. Since the Gerig et al. method did not describe selection 

20 of parameter /c, the following was used: k= V2 cr= ^/l ay/in Equation 10 - the same 

parameter that was used for scale-based diffusion. Except for the result in FIG. 2(d), all results 
for anisotropic diffusion were obtained by applying the method for 3 iterations as 
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recommended in Gerig et a/., 1992, which is herein incorporated by reference. Except for the 
result in FIG. 2(f), all results for scale-based diffusion were obtained by applying the method 
for 7 iterations. 

In Example 1, the three methods are qualitatively compared using a phantom image and 
5 MR images of the heads of multiple sclerosis patients. A quantitative comparison based on a 
phantom containing different geometric shapes is presented in Example 2, to more objectively 
assess the relative performance of each of the three methods. 

Example 1 - Qualitative Comparisons. 
10 In this example, two case studies are presented from MR patient studies, and then two 

phantom experiments are described. 
Comparative Case Study 1.1. 

In FIG. 4, the filtered output from the three methods is illustrated for two regions of 
interest (ROI) - one from the top-left comer and the other from the bottom-right comer - 

15 selected from the MR slice displayed in FIG. 1(a). The two ROIs, and the results from the 

three filtering methods, are shown in FIGs. 4(a)-4(h). FIGs. 4(c)-4(h) depict filtered scenes (in 
2D) corresponding to the scenes in FIGs. 4(a) and 4(b) resulting from anisotropic difftision 
(shown in FIGs. 4(c) and 4(d)), from scale-based neighborhood averaging (shown in FIGs. 4(e) 
and 4(f)), and from scale-based diffiision (shown in FIGs. 4(g) and 4(h)). As can be seen from 

20 all six filtered images (FIGs. 4(c)-4(h)), the smoothness in homogeneous regions seems 

visually to be the same in all images, while boundary contrast seems to be higher using scale- 
based methods, as opposed to using anisotropic diffusion. Also, scale-based methods were 
seen to preserve fine details (for example, the linear bright structiu*e in the left, middle part of 
FIG. 4(b)) better than anisotropic diffiision. Moreover, among the two scale-based methods, 

25 differences did not appear to be visually significant. 
Comparitive Case Study 1 .2. 

FIG. 5 illustrates 3D-comparative results among the three methods. All methods were 
applied on a proton density-weighted 3D MR scene of the head of a multiple sclerosis patient. 
The out-of-plane resolution was seen to be almost 3.5 times more coarse than the in-plane 
30 resolution. The results, thus, also demonstrate the behavior of the three methods under 
resolution anisotropy. 
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The size of the scene domain was 181x246x55, with a voxel size of 0.86x0.86x3.0 
mm^. In the non-optimized experimental implementation, the computation of scale took 10-12 
seconds on a 450 MHz Pentium PC for the scene in FIG. 5, while each iteration of anisotropic 
and scale-based diffusion required 7-8 seconds. Compared to the anisotropic diffusion method, 
5 scale computation was clearly an extra step for the scale-based diffusion method. For the same 
scene, scale-based neighborhood averaging method took 17-18 seconds in total on the same 
machine, which is faster than anisotropic diffusion, which required 21-24 seconds when run for 
three (3) iterations. Moreover, the same improvement of SNR values was achieved in 
homogeneous regions whether anisotropic or scale-based diffusion was applied for the same 

10 number of iterations. However, more blurring was seen across boundaries and in fine 
structured regions in the anisotropic diffusion method. 

Since it was not clear how to best display such 3D scenes (short of showing all slices at 
sufficiently high magnification), one slice is shown in FIG. 5(a) taken from the original scene 
with an ROI. A portion of the sHce was zoomed up (magnified) in FIG. 5(b), to allow a closer 

15 scrutiny of the images. FIGs, 5(c)-(e) show the same ROI for the filtered scenes obtained using 
anisotropic diffusion (FIG. 5(c)), using scale-based neighborhood averaging (FIG. 5(d)), and 
using scale-based diffusion (FIG. 5(e)) methods, respectively. While there are no visually 
significant differences among the three methods in the level of smoothing within homogeneous 
regions, it can been seen that boundary contrast is significantly improved using the scale-based 

20 methods. Some of the lesions (hyper intense blobs) appear blurred in the filtered scene as a 

result of anisotropic diffusion, while by comparison, their contrast is retained using scale-based 
methods. 

Although there is major resolution by anisotropy, the scale-based methods were able to 
retain boundary contrast and subtle structures. As between the two scale-based methods, the 
25 resuh using scale-based diffusion appears visually more pleasing than the result produced using 
scale-based neighborhood averaging. 

Phantom Example 1.1. 

FIG. 2 demonstrates how diffusion proceeds across certain blurred object boundaries in 
the method of anisotropic diffusion (FIGs. 2(c) and 2(d)), whereas it is arrested in the scale- 
30 based method (FIGs. 2(e) and 2(f)). The 2D scene, shown in FIG. 2(b) was a test geometric 
phantom scene created after blurring and adding a zero mean Gaussian noise and a slow 
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varying (ramp) background component across colunms to imitate white matter and gray matter 
regions in the human brain. One sUce was selected from a 3D proton density MR scene of the 
head of a multiple sclerosis patient undergoing routine clinical evaluation. Subsequently, the 
white matter and gray matter regions in the original scene were carefully segmented using a 
5 user-steered segmentation method (Falcao et al. Graphical Models Image Processing 60:233- 
260 (1998)). To the white matter region in the scene, a value was assigned equal to the average 
of the intensities within the white matter region in the original MR slice, and to the background 
region within the brain, a value was assigned equal to the average of the intensities within the 
gray matter region in the original MR slice. 

10 From this starting scene, the 2D phantom scene was generated by applying a Gaussian 

blurring and a zero-mean Gaussian noise, and finally by adding a fixed slow varying (ramp) 
background component across columns. FIGs. 2(a) and 2(b) represent the 2D scenes before 
and after adding noise, blurring, and background variation, respectively. The size of the scene 
was 185x232 with an in-plane resolution of 0.86x0.86 mm^. FIGs. 2(d) and 2(f) were obtained 

15 by applying the anisotropic diffusion and scale-based diffusion methods, respectively. 

In both cases, 10 iterations were used. A few more iterations than the recommended 
number were used for the anisotropic method to emphasize the diffusion across boundaries. On 
the other hand, when the process was run for fewer iterations, /.e., the recommended number of 
3 iterations, the noise reduction was not significant, as can be seen from FIG. 2(c). It can be 

20 observed from FIGs. 2(d) and 2(f), that while the interior (homogeneous) regions were 
(visually) equally smooth in the two images, contrast at the boundaries is significantly 
improved in FIG. 2(f). A similar phenomenon, although to a lesser degree, can be observed 
between FIGs. 2(c) and (e). 

It is unreasonable to expect that a filtering method will retrieve structures that are lost 

25 by any degradation in the imaging process. Therefore, some finer details seen in FIG. 2(a) 
were lost in FIG. 2(f). Given the degraded phantom of FIG. 2(b) as the starting scene, it 
appeared in the filtered scene in FIG. 2(f) that someone had carefully smoothed the scene while 
preserving the structures and boundaries as much as possible. However, any such action would 
require at least a rough knowledge of the object size at various points, which is exactly what is 

30 provided by scale. 
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Phantom Example 1.2. 

The second phantom example was also a 2D phantom containing different geometric 
objects with varying shape, size and orientation. The entire phantom scene domain was divided 
into six square regions of equal size, three of which are at the top, and the other three are at the 
5 bottom. The initial phantom consisted of only two different intensity regions representing the 
objects (brighter) and the background (darker). This was generated by filling each of the six 
square regions with geometric objects of different size, shape and orientation. 

From this initial binary scene, 25 phantom scenes were generated by blurring (using a 
2D Gaussian kernel) the scene at 5 different degrees, and by adding a zero-mean Gaussian- 

10 correlated noise component at 5 levels. Each of the 25 scenes was further modified by adding a 
slow varying (ramp) background component which changed from 0 at the first colunm to 200 at 
the last column. Each of the resulting 25 scenes was filtered using the three filtering methods. 
One of the final phantom scenes, corresponding to a low level of blurring and a medium level 
of noise, is shown in FIG. 6(a). The size of the phantom scene was 800x600. 

15 Due to the large size of the phantom, subtle details are not visible in the display of FIG. 

6(a). The phantom scene, as well as the results using the three filtering methods, are shown in 
FIG. 7 over each of the six ROIs marked in FIG. 6(a). FIG. 7(a) depicts scenes within the six 
ROIs taken from the unfiltered phantom scene, while FIGs. 7(b)-7(d) depict the corresponding 
filtered scenes that resulted fi"om anisotropic diffusion ((shown in FIG. 7(b)), from scale-based 

20 neighborhood averaging (shown in FIG. 7(c)), and fi'om scale-based diffusion (shown in FIG. 
7(d)). 

FIG. 6 also compares the three filtering methods via the intensity profile taken along a 
horizontal line, as shown in FIG. 6(a). FIG. 6(b) shows in magnification the profile over a 
small region indicated at the bottom of FIG. 6(a) (see rectangle outlined in black), and FIGs. 

25 6(c)-(f) show four profiles corresponding to the same region in four different scenes - the initial 
phantom scene after adding only the ramp component (but not adding noise and blurring), and 
the filtered scenes obtained using anisotropic diffusion (shown in FIG. 6(d)), using scale-based 
neighborhood averaging (shown in FIG. 6(e)), and using scale-based diffusion (shown in FIG. 
6(f)), respectively. In FIG. 6(d) (anisotropic filtering), the identities of the five (5) narrow, 

30 vertical colunms, each representing a small square (see FIG. 7, 4th column), were lost; whereas 
by comparison, they were preserved by the scale-based methods (FIGs. 6(e) and 6(f)). 
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Meanwhile, filtering over each piecewise homogeneous region was also improved using 
scale-based diffusion, as compared with anisotropic diffusion. This was also visible in the 
background region shown in the ROI scenes. The first and the sixth columns in FIG. 7 were 
used to demonstrate the behavior of the three methods in homogeneous regions, while other 
5 columns demonstrate their behaviors in regions containing fine structures. Although all 
methods preserve large structures, the scale-based methods produced significantly improved 
results for finer structures. 

Example 2 - Quantitative Comparisons. 
10 The quantitative comparison utilized the scenes depicted in FIGs. 6 and 7. Let Ci = (C, 

fb) denote the original binary phantom scene (before adding artifacts), and \etE= {Ci\ Q^iQ 
fi), 1 < z < 25} denote the set of 25 phantom scenes generated fi-om as described in Example 
1, Phantom Example 1.2. Let let Exi = {Cxi I Cxi = (Cyfxi), 1 < / < 25} denote the set of scenes 
produced by applying the filtering method x e {ad, sa, sd} to the scenes in E, where ad, sa, sd 
15 correspond to anisotropic diffusion, scale-based neighborhood averaging and scale-based 
4=1 diffusion, respectively. Over the scene domain C, 100 ROIs, C"', \<i< 100 were randomly 

J1 selected, each of size 18x18, from homogeneous regions, such that all the pixels in a ROI had 

the same intensity (object or background) in the initial phantom scene ^. Of the 100 ROIs, 60 
were from background regions, and 40 were from object regions. Additionally, 50 ROIs, C^', 1 
£ 20 < / < 50, were selected from regions containing boundaries. Each C^' was selected such that it 
Si contained pixels from both the object and background regions. Let and 5, denote the set of 

C3 pixels in C^' from the object and background regions, respectively. 

Q For jc e {ad, sa, sd} l<i< 25, and 1 <y < 100, "residual noise" RN^^i was defined as 

the standard deviation of pixel intensities fxi{c) over each homogeneous ROI C^^ in a filtered 
25 scene resulting from method x for the phantom scene at the level of blurring and noise given by 
the index 

Tables 1, 2, 3 list in each cell, the mean and the standard deviation of RN{^- , RN^^^. , and 

^idi y respectively, at each degree of blurring and noise for the 100 ROIs. The degrees of 
blurring and noise increase along colunms and rows, respectively, from the upper left comer. 
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Table 1 . The mean (first entry) and standard deviation (second entry) of RN^^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right while blurring increases along columns fi"om top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blur I 


6.036 


8.438 


10.948 


13.457 


16.009 




0.966 


1.278 


1.632 


1.984 


2.341 


Blur 2 


6.076 


8.498 


10.995 


13.506 


16.056 




0.985 


1.288 


1.637 


1.983 


2.338 


Blur 3 


O.U35 


0.5U3 


1 l.UoO 


13.6oi 


10.0U6 




0.973 


1.298 


1.658 


2.015 


2.395 


Blur 4 


6.008 


8.512 


11.115 


13.691 


16.618 




0.965 


1.300 


1.657 


2.018 


2.401 


Blur 5 


5.987 


8.502 


11.111 


13.694 


16.625 




0.963 


1.299 


1.658 


2.021 


2.401 



Table 2 . The mean (first entry) and standard deviation (second entry) of RN^^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
fi-om left to right, while blurring increases along columns firom top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blur I 


4.655 


6.493 


8.813 


11.118 


13.490 




1.350 


1.491 


1.797 


2.157 


2.562 


Blur 2 


4.683 


6.661 


8.938 


11.244 


13.602 




1.339 


1.611 


1.893 


2.226 


2.598 


Blur 3 


4.603 


6.684 


9.282 


11.896 


15.586 




1.214 


1.638 


1.976 


2.354 


2.728 
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Blur 4 


A C OA 


o.oUo 




1 1 Q1 1 


1D.OU4 




1 170 


1 543 


1 932 


2 371 


2 743 


Blur 5 


4.578 


6.795 


9.404 


11.920 


15.613 




1.161 


1.521 


1.919 


2.376 


2.747 



Table 3 . The mean (first entry) and standard deviation (second entry) of RN^^^- values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right, while blurring increases along columns from top to bottom. 





Noise 1 


Noise L 


Noise 3 


Noise 4 


Noise 5 


Blur I 


4.204 


5.920 


7.954 


10.072 


12.334 




0.923 


1.250 


1.631 


2.038 


2.530 


Blur! 


4.223 


5.966 


7.979 


10.114 


12.401 




0.924 


1.269 


1.646 


2.051 


2.529 


Blur 3 


4.238 


5.987 


8.254 


10.778 


14.774 




0.924 


1.276 


1.715 


2.235 


2.855 


Blur 4 


4.239 


6.100 


8.417 


10.785 


14.785 




0.924 


1.298 


1.747 


2.241 


2.854 


Blur 5 


4.240 


6.102 


8.420 


10.787 


14.797 




0.924 


1.298 


1.748 


2.240 


2.857 



From Tables 1, 2, and 3, it can be seen that at every level of blurring and noise, mean 
residual noise is maximal for anisotropic diffusion, and it is minimal for scale-based diffusion. 
Moreover, at each level of blurring and noise, three paired Mests (Rosner, Fundamentals of 
10 Biostatistics , New York, NY, Duxbury Press, 1995, each over 100 pairs (iWA /WA), xi^y.x, 

y € {ad, sa, sd}^ and 1 <y < 100, of data produced by two different filtering methods, were 
conducted under the null hypothesis that there is no statistically significant difference between 
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the two methods. The hypothesis was rejected at a very high level of confidence {p < 10*^) in 
all tests between a scale-based method and anisotropic diffusion. This confidence level was 
less than 0.3x10"^ for the Mest between the two scale-based methods. 

Obviously, a measure of residual noise alone is not sufficient to characterize the 
effectiveness of a filtering method. For this reason, a measure of "relative contrast" was 

defined, RCi^ for jc e {ad, sa, sd), 1 < / < 25, 1 <y < 50 as follows: 



AL/ • — 1 — , 



(Equation 18) 



10 wherein and a^. denote the mean and standard deviation, respectively, of pixel intensities 

fxi{c) over the object region Oj and Mf[ and crff denote similar entities over the background 
region Bj, Notably, RC gives a measure of object-to-background contrast relative to residual 
noise in each region. 

Tables 4, 5, and 6 list in each cell the mean and the standard deviation of RC^^^- , RC^^^ , 

15 and RC{^^ , respectively, at each degree of blurring and noise for the 50 ROIs. Notably, a 
higher value of i?C signifies more accurate filtering. 

Table 4. The mean (first entry) and standard deviation (second entry) of RC{^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right while blurring increases along columns fi"om top to bottom. 

20 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blurl 


7.983 


6.448 


5.185 


4.048 


3.240 




3.400 


2.076 


1.317 


0.772 


0.478 


Blur! 


4.764 


3.691 


2.933 


2.404 


2.059 




2.048 


0.908 


0.307 


0.345 


0.491 


Blur 3 


2.286 


2.002 


1.784 


1.674 


1.609 




0.348 


0.561 


0.651 


0.682 


0.641 
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Blur 4 


1 0*5 1 

1.833 


1.735 


1 /C1 1 

Loll 


1 Anc 
1.475 


1 y1 1 O 
1.416 








0 776 


0 7S1 


0 696 


Blur 5 


1.638 


1.561 


1.468 


1.365 


1.307 




0.755 


0.749 


0.753 


0.758 


0.700 



Table 5. The mean (first entry) and standard deviation (second entry) of RC^^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right while blurring increases along columns from top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise J 


Blur I 


7.755 


6.340 


5.372 


4.627 


4.010 




2.645 


1.455 


0.860 


0.556 


0.412 


Blur 2 


5.176 


4.262 


3.635 


3.140 


2.705 




1.942 


1.138 


0.692 


0.429 


0.230 


Blur 3 


3.311 


2.679 


2.184 


2.002 


1.860 




0.919 


0.289 


0.244 


0.344 


0.369 


Blur 4 


2.390 


2.116 


1.878 


1.675 


1.590 




0.187 


0.266 


0.402 


0.485 


0.462 


Blur 5 


1.881 


1.780 


1.638 


1.501 


1.426 




0.396 


0.434 


0.495 


0.528 


0.502 



Table 6. The mean (first entry) and standard deviation (second entry) of ^C^^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right while blurring increases along columns from top to bottom. 
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Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


BlurX 


8.033 


6.569 


5.563 


4.808 


4.181 




3.044 


1.721 


1.019 


0.662 


0.485 


Blur! 


5.341 


4.386 


3.738 


3.250 


2.847 




2.238 


1.328 


0.809 


0.503 


0.294 


Blur 3 


3.503 


2.923 


2.321 


2.084 


1.913 




1.075 


0.486 


0.246 


0.381 


0.406 


Blur A 


2.690 


2.226 


1.905 


1.673 


1.583 




0.423 


0.237 


0.445 


0.585 


0.530 




1.902 


1.766 


1.612 


1.481 


1.412 


Blur 5 


0.418 


0.492 


0.577 


0.627 


0.566 



From Tables 4 and 6 it can be seen that, at every level of blurring and noise, the mean 
relative contrast for the scale-based filtering method is higher than that for anisotropic 
5 diffusion. At each level of blurring and noise, three paired /-tests each of 50 pairs of {RC\^ , 
RC^yi \x^y,x,y e {ad, sa, sd) , and 1 <7 < 50, of RC data produced by two different filtering 

methods were conducted under the null hypothesis that there is no statistically significant 
difference between the two methods. 

For the /-test between anisotropic diffusion and scale-based diffusion, the hypothesis 

10 was rejected at a very high level of confidence {p < 0.3x10"^) at every level of blurring and 
noise, except for two (2) cases indicated by {Blur 1, Noise 1) and {Blur 1, Noise 2), in which 
two cases the results were not statistically significant. For the /-test between anisotropic 
diffusion and scale-based neighborhood averaging, the hypothesis was rejected at a very high 
level of confidence (p < 0.28x10'^) at every level of blurring and noise, except for three (3) 

15 cases, which were at minimum blurring and the three lowest noise levels, and in which three 
cases the results were not statistically significant. For the /-test between the two scale-based 
methods, scale-based neighborhood averaging and scale-based diffusion, the hypothesis was 
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rejected at a very high level of confidence {p < 0.6x10"^) at every level of blurring and noise, 
except for six (6) cases, which were at maximum blurring and the four highest noise levels and 
at second maximum blurring and the two highest noise levels, and in which six cases the results 
were not statistically significant. 
5 hi conclusion, qualitative comparisons based on geometric phantoms generated under a 

range of conditions of blurring, noise and background variation, as well as patient MR images, 
indicate significant improvements and clear superiority using the scale-based filtering methods 
of the invention over the prior art anisotropic diffiision methods, for preserving boundary 
sharpness and fine details while suppressing noise. As between the two scale-based methods, 

10 results using scale-based diffusion have been found to be superior, except at extremely low 
blurring and noise where the differences were not statistically significant. Both scale-based 
methods embodying the invention use region-constancy based on local homogeneity to effect 
filtering. However, scale-based neighborhood averaging uses it in a hard fashion, while scale- 
based diffusion uses it in a fiizzy manner. Because of this difference in the basic strategy, 

1 5 results using scale-based diffusion were considered to be more visually pleasing than those 
using scale-based averaging. Since the same trend was observed in the quantitative 
experiments, scale-based diffusion was preferred over scale-based neighborhood averaging, 
although the latter is computationally the least expensive among the three methods. 

Each and every patent, patent application and publication that is cited in the foregoing 

20 specification is herein incorporated by reference in its entirety. 

While the foregoing specification has been described with regard to certain preferred 
embodiments, and many details have been set forth for the purpose of illustration, it will be 
apparent to those skilled in the art that the invention may be subject to various modifications 
and additional embodiments, and that certain of the details described herein can be varied 

25 considerably without departing from the spirit and scope of the invention. Such modifications, 
equivalent variations and additional embodiments are also intended to fall within the scope of 
the appended claims. 
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